
### CHANGE THIS TO THE LOCATION OF THE DATA FILES!! ###
dir1 <- "D:\\Research\\Districting\\rep\\"

#######################################
# first, load library and subroutines #
#######################################

# load libraries
library0 <- function(s)
{
	if(s %in% rownames(installed.packages()) == FALSE) {install.packages(s,repos='http://cran.us.r-project.org')}
	require(s,character.only=TRUE)
}
library0("sandwich") # for robust/clustered standard errors
library0("lmtest") # for robust/clustered standard errors
library0("readxl")
library0("fastmatch")
library0("Formula") # for formulas with multiple parts
library0("stringr")

# print and flush
print0 <- function(s1,s2=logical(0))
{
	if(length(s2) == 0)
	{
		print(s1)
		flush.console()
	} else {
		print(paste(s1,": ",s2,sep=""))
		flush.console()
	}
}

# used by sum0 function
apply0 <- function(func1,x,subset=numeric(0),by=numeric(0),bys=numeric(0),...)
{
	if(length(subset) == 0) {
		x1 <- subset(x,!is.na(x))
		if(length(by) != 0) by1 <- subset(by,!is.na(x))
	} else if(length(subset) != length(x)) {
		# Error
		stop("Incompadible Dimensions")
	} else {
		x1 <- subset(x,subset&!is.na(x))
		if(length(by) != 0) by1 <- subset(by,subset&!is.na(x))
	}
	if(length(x1) == 0) {
		NA
	} else {
		if(length(by) == 0) {
			func1(x1,...)
		} else {
			if(length(bys) == 0) bys <- sort(unique(by1))
			dat1 <- cbind(x1,by1)[order(by1),]
			J <- length(unique(by1))
			start <- rep(NA,J)
			start[1] <- 1
			tab1 <- table(by1)
			if(J > 1) {
				for(j in 2:J) start[j] <- start[j-1] + tab1[j-1]
				end <- c(start[2:J]-1,length(x1))
			} else {
				end <- length(x)
			}
			ret <- rep(NA,length(bys))
			for(j in 1:J) {
				k <- fmatch(dat1[start[j],2],bys)
				if(!is.na(k)) ret[k] <- func1(as.numeric(dat1[start[j]:end[j],1]),...)
			}
			names(ret) <- bys
			ret
		}
	}
}

# sum
sum0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	apply0(sum,x,subset,by,bys)
}

# plot points, on a subset
points0 <- function(x,y,subset=logical(0),type='p',col='black',addreg=FALSE,pch=1,lty=1,regcol='black')
{
	if(length(subset) == 0) {
		points(x,y,type=type,col=col,pch=pch,lty=lty)
		if(addreg) { abline(lm(y ~ x),col=regcol) }
	} else {
		x0 <- subset(x,subset)
		y0 <- subset(y,subset)
		points(x0,y0,type=type,col=col,pch=pch,lty=lty)
		if(addreg) { abline(lm(y0 ~ x0),col=regcol) }
	}
}

# plot text, on a subset
text0 <- function(x,y,z,subset=logical(0),col='black',addreg=FALSE,regcol='black',cex=1)
{
	if(length(subset) == 0) {
		text(x,y,z,col=col,cex=cex)
		if(addreg) { abline(lm(y ~ x),col=regcol) }
	} else {
		x0 <- subset(x,subset)
		y0 <- subset(y,subset)
		if(length(z) == length(x)) {
			z0 <- subset(z,subset)
		} else {
			z0 <- z
		}
		if(length(cex) == length(x)) {
			cex0 <- subset(cex,subset)
		} else {
			cex0 <- cex
		}
		text(x0,y0,z0,col=col,cex=cex0)
		if(addreg) { abline(lm(y0 ~ x0),col=regcol) }
	}
}

# compute density, on a subset
density0 <- function(x,bw=logical(0),subset=logical(0),w=logical(0))
{
	if(length(w) == 0) {
		if(length(subset) == 0) {
			x1 <- x[!is.na(x)]
		} else {
			x1 <- x[!is.na(x)&!is.na(subset)&subset]
		}
		if(length(bw) == 0) {
			density(x1)
		} else {
			density(x1,bw)
		}
	} else {
		if(length(subset) == 0) {
			x1 <- x[!is.na(x) & !is.na(w)]
			w1 <- w[!is.na(x) & !is.na(w)]
		} else {
			x1 <- x[!is.na(x)&!is.na(w)&!is.na(subset)&subset]
			w1 <- w[!is.na(x)&!is.na(w)&!is.na(subset)&subset]
		}
		w1 <- w1 / sum(w1)
		if(length(bw) == 0) {
			density(x1,weights=w1)
		} else {
			density(x1,bw,weights=w1)
		}
	}	
}

# define kernels
ker_norm <- function(x) { dnorm(x); } # normal kernel
ker_unif <- function(x) { ifelse(abs(x)<1, 0.5, 0); } # uniform kernel
ker_epa <- function(x) { ifelse(abs(x)<1, (1-x^2)*3/4, 0); } # epanechnikov kernel
ker_tri <- function(x) { ifelse(abs(x)<1, 1-abs(x), 0); } # triangle kernel

set_kernel <- function(ker)
{
	if(ker == "norm") {
		mu2 <- 1
		nu2 <- 1 / (2 * pi^.5)
		kerfunc <- ker_norm
	} else if(ker == "unif") {
		mu2 <- 1 / 3
		nu2 <- 1 / 2
		kerfunc <- ker_unif
	} else if(ker == "tri") {
		mu2 <- 1 / 6
		nu2 <- 2 / 3
		kerfunc <- ker_tri
	} else if(ker == "epa") {
		mu2 <- 1 / 5
		nu2 <- 3 / 5
		kerfunc <- ker_epa
	} else {
		stop("error -- invalid kernel type")
	}
	list(mu2=mu2,nu2=nu2,kerfunc=kerfunc)
}

# kernel regression (locally constant)
kerreg <- function(x,y,ker="norm",subset=logical(0),w=logical(0),h=NA,inference=FALSE,I=201,xgrid=logical(0),cluster=logical(0))
{
	kerstuff <- set_kernel(ker)
	mu2 <- kerstuff$mu2
	nu2 <- kerstuff$nu2
	kerfunc <- kerstuff$kerfunc

	subset[is.na(subset)] <- FALSE
	if(length(w) == 0) {
		if(length(subset) == 0) {
			samp <- !is.na(x) & !is.na(y)
		} else {
			samp <- !is.na(x) & !is.na(y) & subset
		}
	} else {
		if(length(subset) == 0) {
			samp <- !is.na(x) & !is.na(y) & !is.na(w)
		} else {
			samp <- !is.na(x) & !is.na(y) & !is.na(w) & subset
		}
	}
	x1 <- x[samp]
	y1 <- y[samp]
	if(length(w) > 0) w1 <- w[samp]
	if(length(cluster) > 0) cluster1 <- cluster[samp]
	N <- length(x1)
	if(is.na(h)) h <- nrr(x1)
	if(length(xgrid) == 0) xgrid <- min(x1) + (max(x1) - min(x1)) * (0:(I-1)) / I
	ker1num <- matrix(rep(0,N*I),N)
	ker1den <- matrix(rep(0,N*I),N)
	if(length(w) == 0) {
		for(n in 1:N) ker1num[n,1:I] <- y1[n] * kerfunc((xgrid-x1[n])/h)/(N*h)
		for(n in 1:N) ker1den[n,1:I] <- kerfunc((xgrid-x1[n])/h)/(N*h)
	} else {
		for(n in 1:N) ker1num[n,1:I] <- w1[n] * y1[n] * kerfunc((xgrid-x1[n])/h)/(N*h)
		for(n in 1:N) ker1den[n,1:I] <- w1[n] * kerfunc((xgrid-x1[n])/h)/(N*h)
	}
	ker1numsm <- rep(1,N) %*%  ker1num
	ker1densm <- rep(1,N) %*%  ker1den
	ker1est <- ker1numsm / ker1densm
	if(inference)
	{
		# bootstrap 95% CI
		if(length(cluster) == 0) {
			R <- 100 # number of bootstrap replications
			kerestCurr <- matrix(rep(0,R*I),R)
			if(length(w) == 0) {
				XY <- matrix(rep(0,N*2),N)
				XY[1:N,1] <- x1
				XY[1:N,2] <- y1
				for(r in 1:R)
				{
					XYCurr <- XY[sample(nrow(XY),replace=T),]
					ker1num <- matrix(rep(0,N*I),N)
					ker1den <- matrix(rep(0,N*I),N)
					for(n in 1:N) ker1num[n,1:I] <- XYCurr[n,2]*kerfunc((xgrid-XYCurr[n,1])/h)/(N*h)
					for(n in 1:N) ker1den[n,1:I] <- kerfunc((xgrid-XYCurr[n,1])/h)/(N*h)
					ker1numsm <- rep(1,N) %*%  ker1num
					ker1densm <- rep(1,N) %*%  ker1den
					kerestCurr[r,1:I] <- ker1numsm / ker1densm
				}
			} else {
				XYW <- matrix(rep(0,N*3),N)
				XYW[1:N,1] <- x1
				XYW[1:N,2] <- y1
				XYW[1:N,3] <- w1
				for(r in 1:R)
				{
					XYWCurr <- XYW[sample(nrow(XYW),replace=T),]
					ker1num <- matrix(rep(0,N*I),N)
					ker1den <- matrix(rep(0,N*I),N)
					for(n in 1:N) ker1num[n,1:I] <- XYWCurr[n,3]*XYWCurr[n,2]*kerfunc((xgrid-XYWCurr[n,1])/h)/(N*h)
					for(n in 1:N) ker1den[n,1:I] <- XYWCurr[n,3]*kerfunc((xgrid-XYWCurr[n,1])/h)/(N*h)
					ker1numsm <- rep(1,N) %*%  ker1num
					ker1densm <- rep(1,N) %*%  ker1den
					kerestCurr[r,1:I] <- ker1numsm / ker1densm
				}
			}
		} else {
			if(length(w) > 0) stop("error -- this part not done yet!!")

			clusters <- sort(unique(cluster1))
			C <- length(clusters)
			R <- 100 # number of bootstrap replications
			kerestCurr <- matrix(rep(0,R*I),R)
			XY <- list()
			for(c in 1:C) XY[[c]] <- cbind(x1[cluster1==clusters[c]],y1[cluster1==clusters[c]])
			for(r in 1:R)
			{
				csamp <- sample(C,replace=T)
				XCurr <- logical(0)
				YCurr <- logical(0)
				for(c in 1:C) {
					XCurr <- c(XCurr,XY[[csamp[c]]][,1])
					YCurr <- c(YCurr,XY[[csamp[c]]][,2])
				}
				N1 <- length(XCurr)
				ker1num <- matrix(rep(0,N1*I),N1)
				ker1den <- matrix(rep(0,N1*I),N1)
				for(n in 1:N1) ker1num[n,1:I] <- YCurr[n]*kerfunc((xgrid-XCurr[n])/h)/(N1*h)
				for(n in 1:N1) ker1den[n,1:I] <- kerfunc((xgrid-XCurr[n])/h)/(N1*h)
				ker1numsm <- rep(1,N1) %*%  ker1num
				ker1densm <- rep(1,N1) %*%  ker1den
				kerestCurr[r,1:I] <- ker1numsm / ker1densm								
			}
		}
		ylower <- rep(I,0)
		yupper <- rep(I,0)
		for(i in 1:I)
		{
			ylower[i]=quantile(kerestCurr[1:R,i],probs=.025,type=4,na.rm=TRUE)
			yupper[i]=quantile(kerestCurr[1:R,i],probs=.975,type=4,na.rm=TRUE)	
		}
		list(x=xgrid,y=ker1est,ylower=ylower,yupper=yupper)
	} else {
		list(x=xgrid,y=ker1est)
	}
}

# plot confindence intervals with transparency
plotci <- function(dat,col='black',noci=FALSE,lty=1)
{
	if(col == 'black') {
		col2 <- rgb(0,0,0,0.4)
	} else if(col == 'blue') {
		col2 <- rgb(0,0,1,0.4)
	} else if(col == 'purple') {
		col2 <- rgb(1,0,1,0.4)
	} else if(col == 'red') {
		col2 <- rgb(1,0,0,0.4)
	} else if(col == 'orange') {
		col2 <- rgb(1,165/255,0,.4)
	} else if(col == 'magenta') {
		col2 <- rgb(1,0,1,.4)
	} else if(col == 'cyan') {
		col2 <- rgb(0,1,1,.4)
	} else if(col == 'lightblue') {
		col2 <- rgb(0,0,.7,0.2)
		col <- 'blue'
	} else if(col == 'lightpurple') {
		col2 <- rgb(.7,0,.7,0.2)
		col <- 'purple'
	} else if(col == 'lightred') {
		col2 <- rgb(.7,0,0,0.2)
		col <- 'red'
	} else {
		stop("error in plotci -- invalid color")
	}
	points(dat$x,dat$y,col=col,type='l',lty=lty)
	if(!noci & !is.null(dat$ylower) & !is.null(dat$yupper)) polygon(c(rev(dat$x),dat$x),c(rev(dat$ylower),dat$yupper),col=col2,border=NA)
}

# robust standard errors
robust0 <- function(fm)
{
	require(sandwich, quietly = TRUE)
	require(lmtest, quietly = TRUE)
	N <- length(fm$residuals)
	M <- N
	K <- fm$rank
	dfc <- (M/(M-1))*((N-1)/(N-K))
	uj  <- apply(estfun(fm),2, function(x) tapply(x, 1:N, sum));
	vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
	sm <- summary(fm)
	model <- coeftest(fm, vcovCL)
	if(length(sm$r.squared) > 0) {
		list(model=model,N=N,r.squared=sm$r.squared,coef=model[1:K],theta=model[1:K],se=(diag(vcovCL))^.5,V=vcovCL)
	} else {
		list(model=model,N=N,r.squared=fm$r.squared,coef=model[1:K],theta=model[1:K],se=(diag(vcovCL))^.5,V=vcovCL)
	}
}

# determine if a row of a matrix has any missing values
rowmiss <- function(X)
{
	#!apply(X, 1, function(x)!any(is.na(x)))
	rowSums(is.na(X))>0
}

# get the model matrix, including NAs, from a multipart formula
model.matrix0 <- function(f1,rhs=logical(0))
{
	mf <- model.frame(f1,na.action=na.pass)
	if(length(rhs) == 0) mm <- model.matrix(f1,data=mf)
	else                 mm <- model.matrix(f1,data=mf,rhs=rhs)
	mm
}

# get the model response, including NAs, from a multipart formula
model.response0 <- function(f1,lhs=logical(0))
{
	mf <- model.frame(f1,na.action=na.pass)
	if(length(lhs) == 0) mr <- model.response(data=mf)
	else                 mr <- model.part(f1,data=mf,lhs=lhs)
	mr
}

# get all the y's and X's from a multipart formula, eliminating NAs and subseting
formula.data <- function(f,subset=logical(0),list1=logical(0),se="",nointercept=FALSE)
{
	f1 <- Formula(f)
	len1 <- length(f1)
	if(length(len1) != 2) stop("error in formula.data -- invalid formula")
	y <- list()
	X <- list()
	for(i in 1:len1[1]) y[[i]] <- model.response0(f1,lhs=i)
	for(i in 1:len1[2]) X[[i]] <- model.matrix0(f1,rhs=i)
	ptm <- proc.time()
	y0 <- y
	X0 <- X
	n <- nrow(as.matrix(y[[1]]))
	drop <- rep(0,n)
	for(i in 1:len1[1]) drop <- drop | is.na(y[[i]])
	for(i in 1:len1[2]) drop <- drop | rowmiss(X[[i]])
	if(length(subset) == n) drop <- drop | !subset
	else if(length(subset) != 0) stop("error in formula data")
	for(i in 1:len1[1]) y[[i]] <- y[[i]][!drop]
	for(i in 1:len1[2]) X[[i]] <- X[[i]][!drop,]
	if(nointercept) for(i in 1:len1[2]) X[[i]] <- as.matrix(X[[i]][,2:ncol(X[[i]])])
	if(length(list1) > 0) for(i in 1:length(list1)) if(length(list1[[i]]) > 0)
	{
		if(is.matrix(list1[[i]])) {
			if(nrow(list1[[i]]) != n) stop("error in formula.data -- element of list1 does not have n rows")
			list1[[i]] <- as.matrix(list1[[i]][!drop,])
		} else if(is.array(list1[[i]])) {
			if(length(dim(list1[[i]])) != 3) stop("error in formula.data -- array is not a 3D array")
			if(dim(list1[[i]])[1] != n) stop("error in formula.data -- element of list1 does not have n rows")
			list1[[i]] <- as.array(list1[[i]][!drop,,])
		} else {
			if(length(list1[[i]]) != n) stop("error in formula.data -- element of list1 is not of length n")
			list1[[i]] <- list1[[i]][!drop]
		}
	}

	# code to drop clusters for panel Newey-West
	if(se == "pnw" | se == "panelneweywest")
	{
		if(length(list1$unitid) == 0) stop("error in formula.data -- panel newey west standard errors requested, but unitid missing")
		if(length(list1$timeid) == 0) stop("error in formula.data -- panel newey west standard errors requested, but timeid missing")
		N <- length(y[[1]])
		miss <- rep(FALSE,N)
		for(i in 1:len1[1]) miss <- miss | is.na(y[[i]])
		for(i in 1:len1[2]) miss <- miss | is.na(X[[i]])
		unitids <- sort(unique(list1$unitid))
		T <- max(list1$timeid) - min(list1$timeid) + 1
		lag <- floor(4 * (T / 100)^(4/9))	
		pnwuse <- rep(TRUE,N)
		for(i in 1:length(unitids))
		{
			timecurr <- list1$timeid[list1$unitid==unitids[i]]
			misscurr <- miss[list1$unitid==unitids[i]]
			for(j in 1:length(timecurr)) if(misscurr[j]) timecurr[j] <- NA
			use <- TRUE
			for(j in 1:lag) if(sum(!is.na(fmatch(timecurr-j,timecurr))) < 20) use <- FALSE
			if(!use)
			{
				print(paste("Dropping cluster: ",unitids[i],sep=""))
				pnwuse[list1$unitid==unitids[i]] <- FALSE
			}
		}
		for(i in 1:len1[1]) y[[i]] <- y[[i]][pnwuse]
		for(i in 1:len1[2]) X[[i]] <- X[[i]][pnwuse,]
		list1$timeid <- list1$timeid[pnwuse]
		list1$unitid <- list1$unitid[pnwuse]
		if(length(y[[1]]) < 20) stop("error in formula.data -- panel newey west standard errors requested, but less than 20 observations are left")
	}

	list(y=y,X=X,y0=y0,X0=X0,list1=list1)
}

# call the appropriate standard error routine
get.ses <- function(m1,se,unitid=logical(0),timeid=logical(0))
{
	     if(se == "classic"                     ) stop("error in get.ses --- don't call with classic")
	else if(se == "cluster"                     ) res <- cluster0(m1,unitid)
	else if(se == "nw" | se == "neweywest"      ) res <- neweywest0(m1,timeid)
	else if(se == "pnw" | se == "panelneweywest") res <- panelneweywest0(m1,timeid,unitid)
	else if(se == "robust")                       res <- robust0(m1) # the default option
	else                                          stop("invalid se option")
}

# linear regression, with subset, no intercept, and cluster options
lm0 <- function(f1,se="robust",w=logical(0),subset=logical(0),nointercept=FALSE,unitid=logical(0),timeid=logical(0))
{
	# get y and x
	dat <- formula.data(f1,subset,list(unitid=unitid,timeid=timeid,w=w),se=se,nointercept=nointercept)
	if(length(dat$y) != 1) stop("error in lm0 -- invalid number of repsonses")
	if(length(dat$X) != 1) stop("error in lm0 -- invalid RHS")
	y1 <- dat$y[[1]]
	x1 <- dat$X[[1]]
	w1 <- dat$list1$w
	unitid1 <- dat$list1$unitid
	timeid1 <- dat$list1$timeid
	
	# estimate the model
	if(length(w1) == 0) {
		lm1 <- lm(y1 ~ 0 + x1)
	} else {
		lm1 <- lm(y1 ~ 0 + x1,weights=w1)
	}
	if(se=="classic") {
		s1 <- summary(lm1)
		res <- list()
		res$coef <- s1$coefficients[,1]
		res$se <- s1$coefficients[,2]
		res$N <- length(y1)
		res$r.squared <- s1$r.squared
		res$ivs <- names(lm1$coefficients)
	} else {
		res <- get.ses(lm1,se,timeid=timeid1,unitid=unitid1)
		res$ivs <- rownames(res$model)
	}

	if(length(unitid1) == 0) res$Clusters <- 1
	else                     res$Clusters <- length(unique(unitid1))
	
	# model predictions
	x <- dat$X0[[1]]
	if(nointercept) x <- x[,2:ncol(x)]
	pred <- x %*% lm1$coef
	XpX <- t(x1) %*% x1
	res$pred <- pred
	res$XpX <- XpX

	# hack to get correct model fit when intercept is included as a variable
	if(length(w1) == 0) {
		lm2 <- lm(y1 ~ x1)
	} else {
		lm2 <- lm(y1 ~ x1,weights=w1)
	}
	res$theta <- res$coef
	res$r.squared <- summary(lm2)$r.squared
	res$sigma <- summary(lm2)$sigma
	res$AIC <- AIC(lm2)
	res$BIC <- BIC(lm2)
	fac <- (res$N-1)/(res$N-ncol(x))
	res$r.squared.adj <- fac * res$r.squared + 1 - fac
	res$ivs <- substr(res$ivs,3,nchar(res$ivs))
	res$ivs[res$ivs=="tercept)"] <- "Intercept"
	if(se=="classic") {
		res$V <- res$sigma^2 * solve(XpX)
		res$predse <- rep(NA,nrow(x))
		for(n in 1:nrow(x)) res$predse[n] <- res$sigma * (1 + t(x[n,]) %*% solve(XpX) %*% x[n,])^.5
		tcrit <- qt(.975,res$N-ncol(x))
		res$predlow <- res$pred - tcrit * res$predse
		res$predhigh <- res$pred + tcrit * res$predse
	}
	res
}

# ground logit log-likelihood
grouped.logit.ll <- function(beta,x,y,Nj)
{
	prob <- plogis(x %*% beta)
	ll <- log(choose(Nj,y)) + y * log(prob) + (Nj - y) * log(1 - prob)
	-sum(ll)
}

# estimate grouped logit
grouped.logit <- function(f1,subset)
{
	# get y and x
	dat <- formula.data(f1,subset)
	if(length(dat$y) != 2) stop("error in grouped.logit -- invalid number of repsonses")
	if(length(dat$X) != 1) stop("error in grouped.logit -- invalid RHS")
	y  <- dat$y[[1]]
	Nj <- dat$y[[2]]
	x  <- dat$X[[1]]

	# estimate the model
	theta0 <- rep(0,ncol(x))
	opt1 <- optim(theta0,grouped.logit.ll,x=x,y=y,Nj=Nj,method="BFGS",hessian=TRUE,control=list(REPORT=1,trace=6,maxit=1000))
	theta <- opt1$par
	hess <- opt1$hess
	V <- solve(hess)
	se <- sqrt(diag(V))
	list(theta=theta,se=se,V=V,N=length(y),ivs=colnames(x),r.squared=NA,Clusters=NA)
}

# round to k digits
dig <- function(x,k=3)
{
	y <- format(round(x,k),nsmall=k)
	for(i in 1:length(y)) {
		if(is.na(x[i])) y[i] <- ""
	}
	y
}

# round to 3 digits
dig3 <- function(x)
{
	dig(x,3)
}

# coefficient and standard error, nicely formated
as.coef <- function(coef,se,k=3)
{
	if(is.matrix(coef)) {
		res <- coef
		for(i in 1:ncol(coef)) res[,i] <- paste(dig(coef[,i],k),stars(coef[,i],se[,i]),sep="")
		res
	} else {
		paste(dig(coef,k),stars(coef,se),sep="")
	}
}

# standard error, nicely formated
as.se <- function(se,k=3)
{
	if(is.matrix(se)) {
		res <- se
		for(i in 1:ncol(se)) res[,i] <- str_replace_all(paste("(",dig(se[,i],k),")",sep=""),"\\(\\)","")
		res
	} else {
		str_replace_all(paste("(",dig(se,k),")",sep=""),"\\(\\)","")
	}
}

# compute normal pval
pval <- function(est,se)
{
	2*pnorm(-abs(est/se))
}

# add stars to table
stars.scalar <- function(coef,se)
{
	if(is.na(coef)) {
		""
	} else if(is.na(se)) {
		"?"
	} else {
		pval1 <- pval(coef,se)
		if(pval1 <= 0.001) {
			"***"
		} else if(pval1 <= 0.01) {
			"**"
		} else if(pval1 <= 0.05) {
			"*"
		} else if(pval1 <= 0.1) {
			"+"
		} else {
			""
		}
	}
}

# add stars, vector version
stars <- function(coef,se)
{
	if(is.matrix(coef))
	{
		m <- nrow(coef)
		n <- ncol(coef)
		res <- coef
		for(i in 1:m) for(j in 1:n) res[i,j] <- stars.scalar(coef[i,j],se[i,j])
	} else {
		m <- length(coef)
		res <- coef
		for(i in 1:m) res[i] <- stars.scalar(coef[i],se[i])
	}
	res
}

# pad strings in vector so that they are of the same length
pad0 <- function(x)
{
	K <- length(x)
	len <- max(nchar(x))
	for(k in 1:K) if(nchar(x[k]) < len) x[k] <- paste(x[k],paste(rep(" ",len-nchar(x[k])),collapse=""),sep="")
	x
}

# print model output nicely
print.models <- function(ms)
{
	J <- length(ms)
	ivs <- ms[[1]]$ivs
	if(J > 1) for(j in 2:J) ivs <- c(ivs,ms[[j]]$ivs)
	ivs <- unique(ivs)
	v1 <- c("","Variable",ivs,"N","R-Squared")
	v2 <- matrix(rep(NA,J*(4+length(ivs))),ncol=J)
	v3 <- matrix(rep(NA,J*(4+length(ivs))),ncol=J)
	for(j in 1:J)
	{
		coefs <- rep(NA,length(ivs))
		ses <- rep(NA,length(ivs))
		for(k in 1:length(ms[[j]]$ivs))
		{
			coefs[fmatch(ms[[j]]$ivs[k],ivs)] <- ms[[j]]$theta[k]
			ses[fmatch(ms[[j]]$ivs[k],ivs)] <- ms[[j]]$se[k]
		}
		if(!is.null(ms[[j]]$r.squared)) v2[,j] <- c(names(ms)[[j]],"Coef.",as.coef(coefs,ses),ms[[j]]$N,dig3(ms[[j]]$r.squared))
		else                            v2[,j] <- c(names(ms)[[j]],"Coef.",as.coef(coefs,ses),ms[[j]]$N,"")
		v3[,j] <- c("","(s.e.)",as.se(ses),"","")	
		v2[,j] <- pad0(v2[,j])
		v3[,j] <- pad0(v3[,j])
		v2[,j] <- paste(v2[,j],v3[,j],sep=" ")
	}
	v1 <- pad0(v1)
	for(i in 1:length(v1)) print0(paste(c(v1[i],v2[i,]),collapse=" "))
}

#########################
# second, load the data #
#########################

# house election data
data1 <- read_excel(paste(dir1,"lsq_redistricting.xlsx",sep=""),sheet="statecd")
n1 <- nrow(data1)
year1 <- data1$year
state1 <- data1$state
demvotes1 <- as.numeric(data1$dem_votes)
repvotes1 <- as.numeric(data1$rep_votes)
pwin1 <- as.numeric(data1$pwin)
stateyear1 <- paste(state1,year1,sep="_")

# state control of government data
data2 <- read_excel(paste(dir1,"lsq_redistricting.xlsx",sep=""),sheet="stateyear")
n2 <- nrow(data2)
year2 <- data2$year
state2 <- data2$state
key2 <- paste(state2,year2,sep="_")
cont2 <- as.numeric(data2$cont) # state control of government
redist_com2 <- data2$redist_com # redistring controlled by committee?
vra2 <- data2$vra
vra2b <- data2$vrab # alternative coding where partial coverage states (NC,NH) are coded as *not* subject to pre-clearance
redist_year2 <- data2$redist_year # year of last redistricting

# control during last redistricting
cont <- cont2[fmatch(paste(state2,redist_year2,sep="_"),key2)]

# commission during last redistricting?
redist_com <- redist_com2[fmatch(paste(state2,redist_year2,sep="_"),key2)]

# vra during last redistricting?
vra <- vra2[fmatch(paste(state2,redist_year2,sep="_"),key2)]
vrab <- vra2b[fmatch(paste(state2,redist_year2,sep="_"),key2)]

# aggregate to get house vote share and house seat share by state
aggdata1D <- aggregate(demvotes1,by=list(state1,year1),FUN=sum,na.rm=TRUE)
aggdata1R <- aggregate(repvotes1,by=list(state1,year1),FUN=sum,na.rm=TRUE)
aggdata2D <- aggregate(pwin1==1,by=list(state1,year1),FUN=sum,na.rm=TRUE)
aggdata2R <- aggregate(pwin1==0,by=list(state1,year1),FUN=sum,na.rm=TRUE)
aggdata3  <- aggregate(demvotes1, by=list(state1,year1),FUN=length)

key1D <- paste(aggdata1D[,1],aggdata1D[,2],sep="_")
key1R <- paste(aggdata1R[,1],aggdata1R[,2],sep="_")
key2D <- paste(aggdata2D[,1],aggdata2D[,2],sep="_")
key2R <- paste(aggdata2R[,1],aggdata2R[,2],sep="_")
key3  <- paste(aggdata3[,1],aggdata3[,2],sep="_")

demvotes2 <- aggdata1D[,3][fmatch(key2,key1D)]
repvotes2 <- aggdata1R[,3][fmatch(key2,key1R)]
demwin2   <- aggdata2D[,3][fmatch(key2,key2D)]
repwin2   <- aggdata2R[,3][fmatch(key2,key2R)]
numdist2  <- aggdata3[,3][fmatch(key2,key3)]

seatsh <- demwin2 / (demwin2 + repwin2)
totseats <- demwin2 + repwin2
votesh <- demvotes2 / (demvotes2 + repvotes2)
seatsh[is.nan(seatsh)] <- NA
votesh[is.nan(votesh)] <- NA

# divide sample by time
decade1 <- year2 >= 1942 & year2 <= 1960
decade2 <- year2 >= 1962 & year2 <= 1980
decade3 <- year2 >= 1982 & year2 <= 2000
decade4 <- year2 >= 2002 & year2 <= 2018
decsamp0 <- year2 >= 1942 & year2 <= 2018 & numdist2 > 1
decsamp1 <- decade1 &  numdist2 > 1
decsamp2 <- decade2 &  numdist2 > 1
decsamp3 <- decade3 &  numdist2 > 1
decsamp4 <- decade4 &  numdist2 > 1
decstring0 <- "1942 - 2018"
decstring1 <- "1942 - 1960"
decstring2 <- "1962 - 1980"
decstring3 <- "1982 - 2000"
decstring4 <- "2002 - 2018"

sizesamp1 <- decsamp0 & numdist2 >= 2 & numdist2 <= 5
sizesamp2 <- decsamp0 & numdist2 >= 6 & numdist2 <= 10
sizesamp3 <- decsamp0 & numdist2 >= 11 & numdist2 <= 20
sizesamp4 <- decsamp0 & numdist2 >= 21
sizestring1 <- "2 - 5 Cong. Dists."
sizestring2 <- "6 - 10 Cong. Dists."
sizestring3 <- "11 - 20 Cong. Dists."
sizestring4 <- ">21 Cong. Dists."

votetrans <- log(votesh / (1 - votesh))
votetrans[is.infinite(votetrans)] <- NA
demcon <- 1 * (cont == 1)
repcon <- 1 * (cont == -1)
demcon2 <- 1 * (cont2 == 1)
repcon2 <- 1 * (cont2 == -1)
splitcon <- 1 * (cont == 0)
votetrans_dem <- votetrans * demcon
votetrans_rep <- votetrans * repcon

demconnv <- demcon & vra == 0
repconnv <- repcon & vra == 0
demconv <- demcon & vra == 1
repconv <- repcon & vra == 1
votetrans_demnv <- votetrans * demconnv
votetrans_repnv <- votetrans * repconnv
votetrans_vra <- votetrans * vra == 1
votetrans_demv <- votetrans * demconv
votetrans_repv <- votetrans * repconv

demconnc <- demcon & redist_com == 0
repconnc <- repcon & redist_com == 0
demconc <- demcon & redist_com == 1
repconc <- repcon & redist_com == 1
votetrans_demnc <- votetrans * demconnc
votetrans_repnc <- votetrans * repconnc
votetrans_com <- votetrans * redist_com == 1
votetrans_demc <- votetrans * demconc
votetrans_repc <- votetrans * repconc

h <- .2 # set the kernel bandwidth
ker <- "epa"

# vote-seat curve by year
demvotesyr <- sum0(demvotes2,by=year2)
repvotesyr <- sum0(repvotes2,by=year2)
voteshyr <- demvotesyr / (demvotesyr + repvotesyr)

demseatsyr <- sum0(demwin2,by=year2)
repseatsyr <- sum0(repwin2,by=year2)
seatshyr <- demseatsyr / (demseatsyr + repseatsyr)

# urban-rural time series
data3 <- read_excel(paste(dir1,"lsq_redistricting.xlsx",sep=""),sheet="year")
year3 <- data3$year
ur_diff3 <- as.numeric(data3$ur_diff)
ur_diff_lean3 <- as.numeric(data3$ur_diff_lean)

######################
# results start here #
######################

# FIGURE 1 --- urbal-rural plot
plot(year3,ur_diff3,type='l',main="Rural-Urban Difference in Rep. PID Advantage",xlab="Year",ylab="",ylim=c(-.1,.5),xlim=c(1938,2018))
points(year3,ur_diff_lean3,type='l',lty=2)
legend(1940,.5,c("Excluding Leaners","Including Leaners"),lty=c(1,2))

# FIGURE 2 --- vote-seat cure by year
plot(0:1,type='n',xlim=c(.3,.7),ylim=c(.3,.7),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
text0(voteshyr,seatshyr,names(demvotesyr),cex=.5,subset=names(voteshyr)>=1942)
abline(0,1,lty=2)
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(lm(seatshyr ~ voteshyr,subset=names(voteshyr)>=1942),col='magenta')
abline(lm(seatshyr ~ voteshyr,subset=names(voteshyr)>=1992),col='cyan')
legend(.286,.7,c("Linear Fit","Linear Fit (1992-present)"),col=c('magenta','cyan'),lty=c(1,1))

# FIGURE 3 --- vote-seat by state, labeled by controled of state government
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
points0(votesh,seatsh,subset=decsamp0&splitcon)
points0(votesh,seatsh,subset=decsamp0&demcon,col='blue')
points0(votesh,seatsh,subset=decsamp0&repcon,col='red')
legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','black'),lty=c(1,1,1))

# FIGURE 4 --- vote-seat curves by control of government
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))

# FIGURE 5 --- by decade

par(mfrow=c(2,2))

kerD <- kerreg(votesh,seatsh,subset=decsamp1&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp1&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp1&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring1)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp2&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp2&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp2&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring2)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp3&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp3&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp3&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring3)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp4&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp4&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp4&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring4)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# FIGURE 6 --- seperate out VRA states
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&vra==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&vra==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&vra==0,h=h,ker=ker,inference=TRUE)
kerV <- kerreg(votesh,seatsh,subset=decsamp0&vra==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
plotci(kerV,col='black')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control","VRA"),col=c('blue','red','purple','black'),lty=c(1,1,1,1))

# FIGURE 7 --- seperate out VRA states by control
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&vra==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&vra==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&vra==0,h=h,ker=ker,inference=TRUE)
kerVD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&vra==1,h=h,ker=ker,inference=TRUE)
kerVS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&vra==1,h=h,ker=ker,inference=TRUE)
kerVR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&vra==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue',noci=TRUE)
plotci(kerS,col='purple',noci=TRUE)
plotci(kerR,col='red',noci=TRUE)
plotci(kerVD,col='blue',lty=2,noci=TRUE)
plotci(kerVS,col='purple',lty=2,noci=TRUE)
plotci(kerVR,col='red',lty=2,noci=TRUE)
abline(h=.5,lty=2)
abline(v=.5,lty=2)
legend(.0,.97,c("Dem. Control, No VRA","Rep. Control, No VRA","Split Control, No VRA","Dem. Control, VRA","Rep. Control, VRA","Split Control, VRA"),col=c('blue','red','purple','blue','red','purple'),lty=c(1,1,1,2,2,2))

# FIGURE 8 --- seperate out commmisions
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerC <- kerreg(votesh,seatsh,subset=decsamp0&redist_com==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
plotci(kerC,col='orange')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control","Commission"),col=c('blue','red','purple','orange'),lty=c(1,1,1,1))

# FIGURE 9 --- seperate out commmisions by control
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&redist_com==0,h=h,ker=ker,inference=TRUE)
kerCD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&redist_com==1,h=h,ker=ker,inference=TRUE)
kerCS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&redist_com==1,h=h,ker=ker,inference=TRUE)
kerCR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&redist_com==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue',noci=TRUE)
plotci(kerS,col='purple',noci=TRUE)
plotci(kerR,col='red',noci=TRUE)
plotci(kerCD,col='blue',lty=2,noci=TRUE)
plotci(kerCS,col='purple',lty=2,noci=TRUE)
plotci(kerCR,col='red',lty=2,noci=TRUE)
abline(h=.5,lty=2)
abline(v=.5,lty=2)
legend(.0,.97,c("Dem. Cont., No Comm.","Rep. Cont., No Comm.","Split Cont., No Comm.","Dem. Cont., Comm.","Rep. Cont., Comm.","Split Cont., Comm."),col=c('blue','red','purple','blue','red','purple'),lty=c(1,1,1,2,2,2))

# FIGURE 10 --- altenative bandwidths
par(mfrow=c(1,3))
h <- .1
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Vote Share",ylab="Seat Share",main="Bandwidth = 0.1")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
#legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))

h <- .2
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Vote Share",ylab="Seat Share",main="Bandwidth = 0.2")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
#legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))

h <- .4
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Vote Share",ylab="Seat Share",main="Bandwidth = 0.4")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
#legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))
h <- .2

# FIGURE 11 --- control of redistricting over time
demconseats <- rep(NA,39)
splitconseats <- rep(NA,39)
repconseats <- rep(NA,39)
onedistseats <- rep(NA,39)

for(i in 0:38)
{
	yearcurr <- 1942 + 2 * i
	demconseats[i+1] <- sum0(numdist2,subset=year2==yearcurr&demcon&numdist2>1)
	splitconseats[i+1] <- sum0(numdist2,subset=year2==yearcurr&splitcon&numdist2>1)
	repconseats[i+1] <- sum0(numdist2,subset=year2==yearcurr&repcon&numdist2>1)
	onedistseats[i+1] <- sum0(numdist2,subset=year2==yearcurr&numdist2==1)
}

plot(0:1,type='n',xlim=c(1940,2020),ylim=c(0,300),xlab="Year",ylab="Number of Seats")
points0(1942 + 2*0:38,demconseats,col='blue',type='l')
points0(1942 + 2*0:38,splitconseats,col='purple',type='l')
points0(1942 + 2*0:38,repconseats,col='red',type='l')
points0(1942 + 2*0:38,onedistseats,col='orange',type='l')
legend(1942,300,c("Democratic Control","Split Control","Republican Control","One District State"),lty=c(1,1,1,1),col=c('blue','purple','red','orange'))

# FIGURE 12 --- magnitude of redistricting effects
kerD1 <- kerreg(votesh,seatsh,subset=decsamp1&demcon,h=h,ker=ker,inference=TRUE)
kerS1 <- kerreg(votesh,seatsh,subset=decsamp1&splitcon,h=h,ker=ker,inference=TRUE)
kerR1 <- kerreg(votesh,seatsh,subset=decsamp1&repcon,h=h,ker=ker,inference=TRUE)

kerD2 <- kerreg(votesh,seatsh,subset=decsamp2&demcon,h=h,ker=ker,inference=TRUE)
kerS2 <- kerreg(votesh,seatsh,subset=decsamp2&splitcon,h=h,ker=ker,inference=TRUE)
kerR2 <- kerreg(votesh,seatsh,subset=decsamp2&repcon,h=h,ker=ker,inference=TRUE)

kerD3 <- kerreg(votesh,seatsh,subset=decsamp3&demcon,h=h,ker=ker,inference=TRUE)
kerS3 <- kerreg(votesh,seatsh,subset=decsamp3&splitcon,h=h,ker=ker,inference=TRUE)
kerR3 <- kerreg(votesh,seatsh,subset=decsamp3&repcon,h=h,ker=ker,inference=TRUE)

kerD4 <- kerreg(votesh,seatsh,subset=decsamp4&demcon,h=h,ker=ker,inference=TRUE)
kerS4 <- kerreg(votesh,seatsh,subset=decsamp4&splitcon,h=h,ker=ker,inference=TRUE)
kerR4 <- kerreg(votesh,seatsh,subset=decsamp4&repcon,h=h,ker=ker,inference=TRUE)

kerD1$x <- c(0,kerD1$x,1)
kerD1$y <- c(0,kerD1$y,1)
kerS1$x <- c(0,kerS1$x,1)
kerS1$y <- c(0,kerS1$y,1)
kerR1$x <- c(0,kerR1$x,1)
kerR1$y <- c(0,kerR1$y,1)

kerD2$x <- c(0,kerD2$x,1)
kerD2$y <- c(0,kerD2$y,1)
kerS2$x <- c(0,kerS2$x,1)
kerS2$y <- c(0,kerS2$y,1)
kerR2$x <- c(0,kerR2$x,1)
kerR2$y <- c(0,kerR2$y,1)

kerD3$x <- c(0,kerD3$x,1)
kerD3$y <- c(0,kerD3$y,1)
kerS3$x <- c(0,kerS3$x,1)
kerS3$y <- c(0,kerS3$y,1)
kerR3$x <- c(0,kerR3$x,1)
kerR3$y <- c(0,kerR3$y,1)

kerD4$x <- c(0,kerD4$x,1)
kerD4$y <- c(0,kerD4$y,1)
kerS4$x <- c(0,kerS4$x,1)
kerS4$y <- c(0,kerS4$y,1)
kerR4$x <- c(0,kerR4$x,1)
kerR4$y <- c(0,kerR4$y,1)

approxD1 <- approx(kerD1$x,kerD1$y,votesh)$y
approxS1 <- approx(kerS1$x,kerS1$y,votesh)$y
approxR1 <- approx(kerR1$x,kerR1$y,votesh)$y

approxD2 <- approx(kerD2$x,kerD2$y,votesh)$y
approxS2 <- approx(kerS2$x,kerS2$y,votesh)$y
approxR2 <- approx(kerR2$x,kerR2$y,votesh)$y

approxD3 <- approx(kerD3$x,kerD3$y,votesh)$y
approxS3 <- approx(kerS3$x,kerS3$y,votesh)$y
approxR3 <- approx(kerR3$x,kerR3$y,votesh)$y

approxD4 <- approx(kerD4$x,kerD4$y,votesh)$y
approxS4 <- approx(kerS4$x,kerS4$y,votesh)$y
approxR4 <- approx(kerR4$x,kerR4$y,votesh)$y

estD <- rep(NA,n2)
estS <- rep(NA,n2)
estR <- rep(NA,n2)
for(i in 1:n2)
{
	if(!is.na(numdist2[i])) {
		if(numdist2[i] == 1) {
			estD[i] <- votesh[i] > .5
			estS[i] <- votesh[i] > .5
			estR[i] <- votesh[i] > .5
		} else {
			if(decsamp1[i]) {
				estD[i] <- approxD1[i]
				estS[i] <- approxS1[i]
				estR[i] <- approxR1[i]
			} else if(decsamp2[i]) {
				estD[i] <- approxD2[i]
				estS[i] <- approxS2[i]
				estR[i] <- approxR2[i]
			} else if(decsamp3[i]) {
				estD[i] <- approxD3[i]
				estS[i] <- approxS3[i]
				estR[i] <- approxR3[i]
			} else if(decsamp4[i]) {
				estD[i] <- approxD4[i]
				estS[i] <- approxS4[i]
				estR[i] <- approxR4[i]
			}
		}
	}
}

demseats <- rep(NA,39)
expdemseats <- rep(NA,39)
expdemseatsD <- rep(NA,39)
expdemseatsS <- rep(NA,39)
expdemseatsR <- rep(NA,39)

for(i in 0:38)
{
	yearcurr <- 1942 + 2 * i
	demseats[i+1] <- sum0(numdist2*seatsh,year2==yearcurr)
	expdemseats[i+1] <- sum0(numdist2*(demcon*estD+splitcon*estS+repcon*estR),year2==yearcurr)
	expdemseatsD[i+1] <- sum0(numdist2*estD,year2==yearcurr)
	expdemseatsS[i+1] <- sum0(numdist2*estS,year2==yearcurr)
	expdemseatsR[i+1] <- sum0(numdist2*estR,year2==yearcurr)
}

plot(0:1,type='n',xlim=c(1940,2020),ylim=c(175,320),xlab="Year",ylab="Number of Seats")
points0(1942 + 2*0:38,435*voteshyr,type='l',col='magenta')
points0(1942 + 2*0:38,expdemseats,col='black',type='l')
points0(1942 + 2*0:38,expdemseatsS,col='purple',type='l')
abline(h=218,lty=2)
legend(1950,320,c("Expected Democratic Seats","Dem. Vote Share","Expected Democratic Seats under Split Control"),lty=c(1,1,1),col=c('black','magenta','purple'))

# FIGURE 13 --- probability of winning, by decade
par(mfrow=c(2,2))

kerD <- kerreg(votesh,1*demcon,subset=decsamp1,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,1*splitcon,subset=decsamp1,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,1*repcon,subset=decsamp1,h=h,ker=ker,inference=TRUE)
d <- density0(votesh,subset=decsamp1,w=numdist2)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Probability",main=decstring1)
polygon(c(d$x, 1 ),  .15*c(d$y,0 ), col=rgb(0,0,0,.2),border=NA)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(v=.5,lty=2)

kerD <- kerreg(votesh,1*demcon,subset=decsamp2,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,1*splitcon,subset=decsamp2,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,1*repcon,subset=decsamp2,h=h,ker=ker,inference=TRUE)
d <- density0(votesh,subset=decsamp2,w=numdist2)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Probability",main=decstring2)
polygon(c(d$x, 1 ),  .15*c(d$y,0 ), col=rgb(0,0,0,.2),border=NA)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(v=.5,lty=2)

kerD <- kerreg(votesh,1*demcon,subset=decsamp3,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,1*splitcon,subset=decsamp3,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,1*repcon,subset=decsamp3,h=h,ker=ker,inference=TRUE)
d <- density0(votesh,subset=decsamp3,w=numdist2)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Probability",main=decstring3)
polygon(c(d$x, 1 ),  .15*c(d$y,0 ), col=rgb(0,0,0,.2),border=NA)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(v=.5,lty=2)

kerD <- kerreg(votesh,1*demcon,subset=decsamp4,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,1*splitcon,subset=decsamp4,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,1*repcon,subset=decsamp4,h=h,ker=ker,inference=TRUE)
d <- density0(votesh,subset=decsamp4,w=numdist2)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Probability",main=decstring4)
polygon(c(d$x, 1 ),  .15*c(d$y,0 ), col=rgb(0,0,0,.2),border=NA)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(v=.5,lty=2)

# FIGURE 14 --- magnitude of redistricting effects under complete control
plot(0:1,type='n',xlim=c(1940,2020),ylim=c(150,390),xlab="Year",ylab="Number of Seats")
points0(1942 + 2*0:38,435*voteshyr,type='l',col='magenta')
points0(1942 + 2*0:38,expdemseats,col='black',type='l')
points0(1942 + 2*0:38,expdemseatsD,col='blue',type='l')
points0(1942 + 2*0:38,expdemseatsS,col='purple',type='l')
points0(1942 + 2*0:38,expdemseatsR,col='red',type='l')
abline(h=218,lty=2)
legend(1962,390,c("Expected Seats","Democratic Vote Share","Expected Dem. Seats under Dem. Cont.","Expected Dem. Seats under Split Cont.","Expected Dem. Seats under Rep. Cont."),lty=c(1,1,1,1,1),col=c('black','magenta','blue','purple','red'))

# Figure A.1 --- alternative definitions of preclearence

# seperate out VRA states
par(mfrow=c(1,2))

kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&vra==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&vra==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&vra==0,h=h,ker=ker,inference=TRUE)
kerV <- kerreg(votesh,seatsh,subset=decsamp0&vra==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="NC and NH coded as Pre-clearence States")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
plotci(kerV,col='black')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control","VRA"),col=c('blue','red','purple','black'),lty=c(1,1,1,1))

kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon&vrab==0,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon&vrab==0,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon&vrab==0,h=h,ker=ker,inference=TRUE)
kerV <- kerreg(votesh,seatsh,subset=decsamp0&vrab==1,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="NC and NH not coded as Pre-clearence States")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
plotci(kerV,col='black')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# Figure A.2 --- density by regime

plot(density0(votesh,subset=decsamp0&demcon),col='blue',xlim=c(.3,.7),ylim=c(0,6),xlab="Democratic Vote Share",ylab="Density",main="")
points(density0(votesh,subset=decsamp0&splitcon),col='purple',type='l')
points(density0(votesh,subset=decsamp0&repcon),col='red',type='l')

# Figure A.3 --- density w/ data
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(-.4,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))
points(density0(votesh,subset=decsamp0&demcon)$x,-.35+(1/14)*density0(votesh,subset=decsamp0&demcon)$y,col='blue',xlim=c(.3,.7),ylim=c(0,6),xlab="Democratic Vote Share",ylab="Density",main="",type='l',lty=2)
points(density0(votesh,subset=decsamp0&splitcon)$x,-.35+(1/14)*density0(votesh,subset=decsamp0&splitcon)$y,col='purple',type='l',lty=2)
points(density0(votesh,subset=decsamp0&repcon)$x,-.35+(1/14)*density0(votesh,subset=decsamp0&repcon)$y,col='red',type='l',lty=2)

# Figure A.4 --- density, by decade
par(mfrow=c(2,2))

plot(density0(votesh,subset=decsamp1&demcon)$x,density0(votesh,subset=decsamp1&demcon)$y,col='blue',ylim=c(0,7),main=decstring1,type='l',xlab="Democratic Vote Share",ylab="Density",xlim=c(.3,.7))
points(density0(votesh,subset=decsamp1&splitcon),col='purple',type='l')
points(density0(votesh,subset=decsamp1&repcon),col='red',type='l')

plot(density0(votesh,subset=decsamp2&demcon)$x,density0(votesh,subset=decsamp2&demcon)$y,col='blue',ylim=c(0,7),main=decstring2,type='l',xlab="Democratic Vote Share",ylab="Density",xlim=c(.3,.7))
points(density0(votesh,subset=decsamp2&splitcon),col='purple',type='l')
points(density0(votesh,subset=decsamp2&repcon),col='red',type='l')

plot(density0(votesh,subset=decsamp3&demcon)$x,density0(votesh,subset=decsamp3&demcon)$y,col='blue',ylim=c(0,7),main=decstring3,type='l',xlab="Democratic Vote Share",ylab="Density",xlim=c(.3,.7))
points(density0(votesh,subset=decsamp3&splitcon),col='purple',type='l')
points(density0(votesh,subset=decsamp3&repcon),col='red',type='l')

plot(density0(votesh,subset=decsamp4&demcon)$x,density0(votesh,subset=decsamp4&demcon)$y,col='blue',ylim=c(0,7),main=decstring4,type='l',xlab="Democratic Vote Share",ylab="Density",xlim=c(.3,.7))
points(density0(votesh,subset=decsamp4&splitcon),col='purple',type='l')
points(density0(votesh,subset=decsamp4&repcon),col='red',type='l')

# Figre A.5 --- pre-72 vs. post 72
par(mfrow=c(1,2))

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2<1972&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2<1972&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2<1972&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="Pre-1972")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="Post-1972")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

## Figure A.6 --- 10-year periods

par(mfrow=c(2,4))

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1942&year2<1952&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1942&year2<1952&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1942&year2<1952&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1942-1950")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1952&year2<1962&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1952&year2<1962&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1952&year2<1962&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1952-1960")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1962&year2<1972&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1962&year2<1972&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1962&year2<1972&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1962-1970")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&year2<1982&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&year2<1982&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1972&year2<1982&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1972-1980")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1982&year2<1992&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1982&year2<1992&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1982&year2<1992&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1982-1990")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1992&year2<2002&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1992&year2<2002&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=1992&year2<2002&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="1992-2000")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2002&year2<2012&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2002&year2<2012&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2002&year2<2012&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="2002-2010")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2012&year2<=2022&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2012&year2<=2022&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&year2>=2012&year2<=2022&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="2012-2018")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# Figure A.7 --- by state size
par(mfrow=c(2,2))

kerD <- kerreg(votesh,seatsh,subset=sizesamp1&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=sizesamp1&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=sizesamp1&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=sizestring1)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=sizesamp2&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=sizesamp2&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=sizesamp2&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=sizestring2)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=sizesamp3&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=sizesamp3&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=sizesamp3&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=sizestring3)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=sizesamp4&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=sizesamp4&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=sizesamp4&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=sizestring4)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# Figure A.8 --- unweighted vs. weighted
par(mfrow=c(1,2))

# vote-seat curves by control of government
kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="Unweighted")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))

# weighted, vote-seat curves by control of government
kerD <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main="Weighted")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)
legend(.0,.97,c("Democratic Control","Republican Control","Split Control"),col=c('blue','red','purple'),lty=c(1,1,1))

# Figure A.9 --- by decade, weighted
par(mfrow=c(2,2))

kerD <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp1&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp1&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp1&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring1)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp2&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp2&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp2&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring2)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp3&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp3&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp3&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring3)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp4&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp4&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,w=numdist2,subset=decsamp4&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Democratic Vote Share",ylab="Democratic Seat Share",main=decstring4)
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# Figure A.10 --- clustered ses
par(mfrow=c(1,2))

kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Vote Share",ylab="Seat Share",main="Regular Confidence Bands")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

kerD <- kerreg(votesh,seatsh,subset=decsamp0&demcon,h=h,ker=ker,inference=TRUE,cluster=state2)
kerS <- kerreg(votesh,seatsh,subset=decsamp0&splitcon,h=h,ker=ker,inference=TRUE,cluster=state2)
kerR <- kerreg(votesh,seatsh,subset=decsamp0&repcon,h=h,ker=ker,inference=TRUE,cluster=state2)
plot(0:1,type='n',xlim=c(0,1),ylim=c(0,1),xlab="Vote Share",ylab="Seat Share",main="Clustered Confidence Bands")
plotci(kerD,col='blue')
plotci(kerS,col='purple')
plotci(kerR,col='red')
abline(h=.5,lty=2)
abline(v=.5,lty=2)
abline(0,1,lty=2)

# set up for linear/logit models

demcon_numdist2 <- demcon * numdist2
repcon_numdist2 <- repcon * numdist2

year2adj <- year2  - 1942

demcon_year2 <- demcon * (year2 - 1942)
repcon_year2 <- repcon * (year2 - 1942)

demcon_redist_com <- demcon * redist_com
repcon_redist_com <- repcon * redist_com

demcon_vra <- demcon * vra
repcon_vra <- repcon * vra

year_votetrans <- (year2 - 1942) * votetrans
demcon_year_votetrans <- demcon * (year2 - 1942) * votetrans
repcon_year_votetrans <- repcon * (year2 - 1942) * votetrans

redist_com_votetrans <- redist_com * votetrans
demcon_redist_com_votetrans <- demcon * redist_com * votetrans
repcon_redist_com_votetrans <- repcon * redist_com * votetrans

vra_votetrans <- vra * votetrans
demcon_vra_votetrans <- demcon * vra * votetrans
repcon_vra_votetrans <- repcon * vra * votetrans

# TABLE A.1 --- linear model
lm1 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp0)
lm2 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp1)
lm3 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp2)
lm4 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp3)
lm5 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp4)
print.models(list("All Years"=lm1,"1942-1960"=lm2,"1962-1980"=lm3,"1982-2000"=lm4,"2002-2018"=lm5))

# TABLE A.2 --- linear model, moderators
lm6 <- lm0(seatsh ~ demcon + repcon + year2adj +  demcon_year2 + repcon_year2 + redist_com + demcon_redist_com + repcon_redist_com + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + year_votetrans +  demcon_year_votetrans + repcon_year_votetrans + redist_com_votetrans + demcon_redist_com_votetrans + repcon_redist_com_votetrans + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp0)
lm7 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp1)
lm8 <- lm0(seatsh ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp2)
lm9 <- lm0(seatsh ~ demcon + repcon + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp3)
lm10 <- lm0(seatsh ~ demcon + repcon + redist_com + demcon_redist_com + repcon_redist_com + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + redist_com_votetrans + demcon_redist_com_votetrans + repcon_redist_com_votetrans + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp4)
print.models(list("All Years"=lm6,"1942-1960"=lm7,"1962-1980"=lm8,"1982-2000"=lm9,"2002-2018"=lm10))

# TABLE A.3 --- grouped logit
logit1 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp0)
logit2 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp1)
logit3 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp2)
logit4 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp3)
logit5 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp4)
print.models(list("All Years"=logit1,"1942-1960"=logit2,"1962-1980"=logit3,"1982-2000"=logit4,"2002-2018"=logit5))

# Table A.4 --- grouped logit, moderators
logit6 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + year2adj +  demcon_year2 + repcon_year2 + redist_com + demcon_redist_com + repcon_redist_com + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + year_votetrans +  demcon_year_votetrans + repcon_year_votetrans + redist_com_votetrans + demcon_redist_com_votetrans + repcon_redist_com_votetrans + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp0)
logit7 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp1)
logit8 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + votetrans + votetrans_dem + votetrans_rep,subset=decsamp2)
logit9 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp3)
logit10 <- grouped.logit(demwin2 | totseats ~ demcon + repcon + redist_com + demcon_redist_com + repcon_redist_com + vra + demcon_vra + repcon_vra + votetrans + votetrans_dem + votetrans_rep + redist_com_votetrans + demcon_redist_com_votetrans + repcon_redist_com_votetrans + vra_votetrans + demcon_vra_votetrans + repcon_vra_votetrans,subset=decsamp4)
print.models(list("All Years"=logit6,"1942-1960"=logit7,"1962-1980"=logit8,"1982-2000"=logit9,"2002-2018"=logit10))

